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Abstract 

For systems in nonequilibrium steady states, a novel modulated Gaussian probability distribution 
is derived to incorporate a new phenomenon of biased current fluctuations, discovered by recent 
laboratory experiments and confirmed by molecular dynamics simulations. Our results consistently 
extend Onsager-Machlup fluctuation theory for systems in thermal equilibrium. Connections with 
the principles of Statistical Mechanics due to Boltzmann and Gibbs are discussed. At last, the 
modulated Gaussian distribution is of potential interest for other statistical disciplines, which make 
use of the Large Deviation theory. 
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I. INTRODUCTION 


In this letter the fluctuation theory of currents, originally proposed by Onsager and 
Machlup for systems in thermal equilibrium Ref. [S], is extended to the nonequilibrium 
Steady States (SS). We suggest a novel Modulated Gaussian (MG) probability distribution, 
which supersedes the normal distribution, used to describe the fluctuations of thermody¬ 
namic variables in equilibrium. This new model features an asymmetry of the current 
fluctuations, a property inherent to the nonequilibrium regime. MG includes the Gaussian 
distribution as a special case so that the earlier theory is consistent with ours. Furthermore, 
in Sec. |V| connections between our formalism and the principles of statistical mechanics due 
to Boltzmann and Gibbs are discussed. 

Our research was motivated by the evidence of non-Gaussianity of current fluctuations 
in a nonequilibrium SS, found in recent laboratory experiments Ref. [3] on a shear flow 
in a 2-Dimensional (2D) dusty plasma, similar to that considered in Ref. |5]. The new 
probability distribution, which is based on the Large Deviation (LD) approach (Ref. [ 12 ]), 
allows to account for this newly observed phenomenon. Our theoretical results are confirmed 
by Molecular Dynamics (MD) simulations of a shear flow for a Debye-Hiickel (DH) plasma. 

The asymmetry of nonequilibrium SS current fluctuations, studied in this paper, is due 
to a probability bias, namely due to the higher probability for the current to deviate from 
its most likely value towards larger magnitudes, rather than to smaller ones. Naturally, 
equilibrium systems are free of such a bias, since their most likely magnitude of a current is 
zero B 

Finally, we would like to remark, that the MG probability distribution, which is derived in 
Sec. |III| without any specific physical assumptions, relies solely on the LD theory. Therefore, 
whenever the latter holds, the MG should also be applicable, in principle, to random variables 
outside the held of Statistical Mechanics. 


^ The temporal asymmetry predicted for trajectories of fluctuations onset and decay by the so called 
macroscopic fluctuation theory Ref. [2] is different from the bias of the time-independent nonequilibrium 
SS probabilities, studied here. 
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II. COMPUTATIONAL MODEL 


In our numerical experiments we studied an iso-thermal MD model of a 2D plasma system. 
We will concentrate on computations performed with N = 100 particles of a unit mass m, 
interacting through the DH potential ^dh, shifted and truncated at a cutoff Vc- 

eA [r-Wxp(-^)-r-Wxp(-5f)] , if r < r, 

^DH{r) = < , ( 1 ) 

I 0 , if r > Tc 

where r is the distance between two interacting particles, while the energy constant e and 
the Debye screening length A are chosen as the basis of a reduced units system. 

A constant shear rate and temperature of the plasma were maintained by the SLLOD 
equations of motion coupled to the NosAHoover (NH) thermostat in D = 2 dimensions 
(Ref. i): 


Pi , Y 

qi = -h JQiyX 

m 

Pi = Fi{qi) - 'ypiyX - a^HPi 


N 


^nh = 0 ^ i ^ 


Pi 


, 2 = 1 


ONksTm 


- 1 


( 2 ) 


Here X is a unit vector along the Cartesian x-coordinate axis, ks is the Boltzmann 
constant, 7 = ^ is the applied shear rate, i.e. the gradient of the streaming velocity 
u = (Ma; 0 ); while qi, Pi = m{qi — u) and Fi are, respectively, the i-th particle position {qiy 
being its y-coordinate), peculiar momentum and force on the Ath particle, resulting from 
the interactions with all the other particles; finally, the NH thermostat is characterized by 
its temperature T, the relaxation time 6 and the time-dependent coupling a^n- The flux of 
x-momentum along the y-coordinate generated by Eq. (|^ in the primary simulation cell 
of the linear dimension L, is 


N 

i=l 


PixPiy 

m 


FixVi 


( 3 ) 


^ I.e. the xy-component of the pressure tensor 
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In all our simulations is set equal to le, whereas 6 was adopted as the unit of time. 
The equations of motion were integrated by the classical 4th-order Runge-Kutta algorithm 
(Ref. [7]), with the time step 0.0016* and the potential cutoff = 3A. The linear size of 
primary cell L was 8.86A. 


III. MODULATED GAUSSIAN DISTRIBUTION 

In this section we present a derivation of the MG distribution for the flux Eq. ([^. We 
begin with a probability density function of the form: 


p(P,,) = (4) 

Kb 

The function S{Pxy), which depends implicitly on N, will be interpreted in Sec. Con¬ 
sidering each term of the summation operator in Eq. ([^ as a random variable, we assume, 
that their spatial average P^y obeys the LD theory. That is, there exists a rate function 

I{P,y) = - lim {S{P,y)/NkB}, (5) 

N^oo 

which has a global minimum 

HPxy) = 0 ( 6 ) 

at the most likely value of P^y = Pxy [cf. Ref. [ 12 ] )■ 

Expanding S{Pxy) in Eq. (|^ about Pxy in a power series up to a prescribed hnite order 
n, we obtain: 


S{Pxy) 

ks 


^ 5(0 . 

~ kB I'.kB ^ 

1=1 

de^F + ^S{^Pxy) 

ks 



( 7 ) 


which dehnes (=^) a constant S and a cost function AS{APxy) of the deviation APxy = 
Pxy — Pxy, which will be considered from a physical perspective in Sec. 

From now on, the values of the derivatives S^^\Pxy) will be denoted by Si, for brevity. 
Since the cost function is zero for Pxy = Pxy, Eqs. (j^ suggest to pose: 
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( 8 ) 


As p{Pxy) is the global maximum of probability density, it follows that the derivative Si 
vanishes. For a hnite N Eq. Q then becomes: 


p{P^y) ^ exp{-— + 
kf 


S , A^(AP,^) 


} = exp{-^ + ^^(AP,j^)*}, 

ks ks ^ i\kB 

1 = 1 


( 9 ) 


where the normalization of the total probability requires that 


S 

exp{ —} = / dP^yAS{Pxy - Pxy)- 

J —oo 

For n = 2, Eq. 0 turns into a Gaussian distribution. One may recognize that, this case 
was treated by Onsager and Machlup in Ref. P Uni for the equilibrium state (Pj,^ = 0 ) 
Therefore to account for non-Gaussian phenomena, one must consider higher order terms of 
the cost function, which can be arranged as follows: 

^ i=3 ' ^ 

where the term in the curly braces = 1 + 2 Yli =3 'i^^^xy'^ ^ modulating factor, to 

which the MG distribution owes its name. 

The original Gaussian is recovered, when = 1 . For p{Pxy) to be integrable in Eq. (|^, 
the order n has to be restricted to even integers. Using the lowest order non-Gaussian 
approximation, i.e. n = 4, we replace the derivatives Si by three parameters of scale If, of 
asymmetry A, and of non-Gaussianity P, which allow a more convenient interpretation, in 
Eq. g: 

^ Especially, current fluctuations about equilibrium are considered in Ref. 0. 
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p{Pxy) oc exp{AS{APxy)} = 


2k 


B 


-T,4^(APxy)} 


= exp{-i^^(l - 2,/^AB^^ + 


( 11 ) 


2n2 ' ' n n2 

The dimensionless non-negative constant B > 0 controls the level of non-Gaussianity, 
thus p{Pxy) is exactly Gaussian for B = 0 with the scale parameter If. A determines 
the asymmetry of p{Pxy), which is skew to the left (right), i.e. has a negative (positive) 


skewness, when A < 0 (A > 0), respectively. The numerical factor 2'\/2/3 in Eq. (11), as 
a coefficient of the term with A, was chosen to make AS{APxy) a non-concave function of 
P^y, i.e. > 0, when —1 < v4 < 1. Violation of the non-concavity condition would 

admit special artifacts of the probability distribution, e.g. a local maximum of its density. 
For experimental or numerical observations, such flexibility may result in a data over-fitting, 
and hence should be suppressed by using constrained optimization techniques. 


IV. NON-GAUSSIANITY OF CURRENT FLUCTUATIONS 

Fig. a illustrates probability density estimates of the momentum current Pxy, at low 
and high shear rates, for our computational experiments, including the fits of the Gaussian 
and modulated Gaussian distributions (Eq. (11)) For the low value of 7 = 0.20“^ an 
excellent agreement is observed between the shape of histograms and the parametric models 
(indistinguishable on the graph for the lower shear rate). However, the Gaussian model 
renders an unsatisfactory result in the case of a larger value 7 = 1.06^“^, where the superior 
performance of the MG distribution is evident. 

Most obviously, the Gaussian model fails, when the fluctuations are asymmetric about 
the mode of the probability distribution Pxy, as can be quantified by the skewness statis¬ 
tics. Increasing the departure from equilibrium augments the non-Gaussianity of p{Pxy), 
as confirmed by the plots of the skewness and excess kurtosis, which are both zero for the 
normal distribution, vs. the magnitude of the shear rate in Fig. This occurs, because in 
^ Everywhere we adopt the Freedman-Diaconis rule (Ref. [5]) to calculate the bin width of traditional 

histograms and the band width of smooth histograms with the Gaussian kernel (Ref. [11)1. 

® Because the normalizing constant of the MG can be evaluated only numerically, to obtain the corre¬ 
sponding Maximum Likelihood Estimator (MLE) of the MG model, we used the Wolfram Mathematica® 
implementation of the interior point method with the constraint of non-concavity, mentioned at the end 
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□ Traditional histogram ■ Smooth histogram 
■ Gaussian model ■ Modulated Gaussian 


FIG. 1. Probability density of Pxy'- (a) for a low shear rate in the left panel; (b) for a high shear 
rate in the right panel. For the low shear rate the Gaussian and the MG are indistinguishable. 

nonequilibrium the probability of a positive deviation AP^y > 0 from the most likely value 
Pxy becomes different from the probability of the opposite deviation —APxy. Fig. also 
demonstrates the efficiency of the modulated Gaussian to £t high-order statistics. 





X Observations of Pxy 
• Modulated Gaussian estimation 


FIG. 2. Statistics of skewness and excess kurtosis as a function of the shear rate for our observa¬ 
tions of Pxy and the corresponding fits by the MG. 

For a positive shear rate, which maintains a negative flux Pxy {cf. Fig. [^, we observe 
p{APxy > 0) < p{—APxy). The MG distribution accounts for this bias of probabilities and 
fits very well the 3rd and 4th order statistics (skewness and kurtosis, respectively). We 
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conjecture, that the asymmetry of the distribution p{Pxy) is a characteristic property of the 
current fluctuations about a nonequilibrium SS. It is small, though, in the near-equilibrium 
SS regime and vanishes only in the limit 7 —)■ 0. 

V. CONNECTION WITH CLASSICAL STATISTICAL MECHANICS 


The results obtained in Sec. |III[ can be readily connected with the Boltzmann and Gibbs 


principles of statistical mechanics. According to the Boltzmann principle, given a measure 
of system microstates w{Pxy), for a given value of Pxy, and the total measure of the ac¬ 
cessible microstates W = w{Pxy)dPxy under the specihed macroscopic constraints of 
temperature, volume, shear rate etc., the SS probability density p{Pxy) and the Boltzmann 
entropy SsiPxy) are given by: 



Ssi.Pxy') In la (fAj/)) 


( 12 ) 


respectively. 


Using the notion of total entropy Stot = fcslnlU after Ref. [T], we deduce from Eq. ( 12 ) 


that: 


ks lnp{Pxy) = ks \n.w{Pxy) - /cb InlU = SsiPxy) - SsiPxy) + SsiPxy) - Stot 
= ASsiPxy) + SsiPxy) - Stot, 


(13) 


where in the second equality we added and subtracted Ss^Pxy) to introduce the Boltzmann 
entropy difference ASsiPxy) = Sb^Pxv) - SB^Pxy)- 


Comparing Eq. 0 with Eq. ([I^, one sees that 



( 14 ) 


because ASb^Pxv) AS{Px;y — Pxy) are both zero at P^y = Pxy by dehnition. 





Hence Eq. (14) provides the interpretation of the cost fnnction AS{APxy) and the con¬ 
stant S, introduced in Sec. Ill, in terms of Boltzmann entropy and the total entropy. The 


most likely value of P^y = Pxy maximizes the Boltzmann entropy (c/. Eq. (12)). Consis¬ 
tently, AS{APxy) < 0 is the (Boltzmann) entropy cost of a fluctuation P^y = Pxy + APxy 
1^ The constant S is the remaining total entropy, after subtracting the entropy cost of the 
most likely macrostate. 

Finally, the Gibbs entropy, given by a functional 5 'g[p], for the distribution Eq. 0 is 


>S'g[p(-)] = J^^dPxyp{Pxy)lrip{Pxy) = -kB{lnp{Pxy)) 

= {S - AS (Pxy)) = S - {AS (Pxy)) 

= Stot - SeiPxy) - {^S{Pxy)), (15) 

which extracts 5, up to a constant term {AS{Pxy)), from the SS probability density 
p{Pxy)- 


VI. CONCLUSION 


In Sec. Ill we used the 4th order modulating factor (S 4 ) to derive our modulated Gaussian 
probability distribution. In principle, the same procedure can be applied to obtain the 
next order approximation, using Eg. However, the formalism becomes then much more 


complicated. The performance of the forth order factor was demonstrated in Sec. [IV] and 
proved to be sufficient to describe the asymmetry of ffuctuations and the decay of their 
probability distribution tails. Although we presented our data only for the simulations with 
V = 100 particles, we checked that the modulated Gaussian performs equally well for smaller 
and larger values N as well. 

An important consequence of the asymmetric fluctuations about the SS is the discrepancy 
between the current average and its most likely value {Pxy) 7 ^ Pxy, which are equal in equi¬ 
librium. Since we do not provide the dynamics of fluctuations, like the Langevin equation 
suggested by Onsager and Machlup for equilibrium, it is an open question, how this prob¬ 
ability bias emerges. There are various ways to modify the Langevin equation to produce 


One can also define the (free) energy cost of a fluctuation AF = AS{Pxy)T, cf. Ref. [S]. 
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asymmetric time-independent distribution, e.g. non-Gaussian noise or non-linear determin¬ 
istic term. Physically, this could be a consequence of the arrow of time in a nonequilibrium 
stystem, which admits a bias. Yet a dynamical theory remains to be found. 
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